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I present a fast algorithm to find apparent horizons. This algorithm uses an explicit representation 
of the horizon surface, allowing for arbitrary horizon resolutions and, in principle, shapes. Novel in 
this approach is that the tensor quantities describing the horizon live directly on the horizon surface, 
yet are represented using Cartesian coordinate components. This eliminates coordinate singularities, 
and leads to an efficient implementation. The apparent horizon equation is then solved as a nonlinear 
elliptic equation with standard methods. I explain in detail the coordinate systems used to store 
and represent the tensor components of the intermediate quantities, and describe the grid boundary 
conditions and the treatment of the polar coordinate singularities. Last I give as examples apparent 
horizons for single and multiple black hole configurations. 



I. INTRODUCTION 

Apparent horizons serve several purposes during the 
numerical evolution of general relativistic spacetimes. In 
vacuum, spacetime has no other prominent features that 
can be determined locally in time, such as the shock 
fronts found in hydrodynamics. Event horizons are global 
structures, and it is not possible to find them during a 
time evolution, because the future of the spacetime is not 
yet known. Yet each apparent horizon indicates, under 
reasonable assumptions, that there is an event horizon 
present somewhere at or outside the apparent horizon. 
Furthermore, apparent horizons are under certain cir- 
cumstances also isolated horizons, making it then pos- 
sible to calculate their mass and spin. 

Locating apparent horizons is also necessary when one 
wants to apply excision boundary conditions in a numeri- 
cal evolution, where the apparent horizon has to be used 
to determine the location and shape of the excised re- 
gions. It is therefore important to have a fast and robust 
apparent horizon finder that is closely integrated with 
the evolution code. 



II. METHOD 

An apparent horizon is an outermost 2-surface within 
a spacelike hypersurface (time slice) with zero expansion 
and spherical topology. The condition that it be an out- 
ermost such surface is difficult to verify in practice, and 
I will not require it in the following. That means that in 
each point of the horizon, the condition 

// : V,,' • A,,: s'.s' .,/'•' , II (1) 

has to hold (see e.g. [15, chapter 4], or [16]). Here s 1 is the 
(spacelike) outward normal to the horizon, and gij and 
Kij are the three-metric and extrinsic curvature of the 
time slice. H is the so-called apparent horizon function. 
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An apparent horizon surface can be represented ex- 
plicitely by a function h(9, </>) that specifies (this is an 
arbitrary choice) the radius r as a function of the spheri- 
cal coordinates 8 and <j>. This choice restricts the possible 
surfaces to those of star-shaped regions about the origin, 
but this restriction does not cause problems in practice; 
other choices are equally possible. The horizon consists 
out of those points (r, 0, (f>) of the time slice that satisfy 
the condition r = h(9, </>). 

The apparent horizon can also be defined implicitly 
through an expression F{x l ) = 0, with the level set func- 
tion F defined e.g. as 

F{r,6,<j>)=r-h{e,<j>) . (2) 

This representation leads to the equation 
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for the outward normal. The function F and the space- 
like normal s z can be calculated using either spherical 
coordinates (r, 9, </>), or transformed to use Cartesian co- 
ordinates [x, y, z). Such coordinate distinctions are actu- 
ally rather important for the numerical implementation 
of an algorithm, so that I do want to treat them in some 
detail. 

The condition H = is a nonlinear elliptic equation in 
h, as can be seen by substituting the definitions of s l and 
F into H. With the above coordinate choice, the domain 
of the equation H(9,(f>) = is (0,0) E (0,tt) x [0, 2tt) 
with coordinate singularities at the poles; the range of 
the equation is h = r € (0, oo). 

As alternatives to treating the apparent horizon equa- 
tion as an elliptic equation, as I do, there are also other 
methods for solving it to be found in the literature. Com- 
monly used methods are mean curvature flow [13] or level 
flow methods [14], or minimisation procedures [2]. The 
advantages of flow methods are that they are in gen- 
eral more robust in finding an horizon, and that they are 
in practice able to find the outermost apparent horizon, 
given a large sphere as initial data. Their disadvantage 
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is that they lead to parabolic or degenerate elliptic equa- 
tions that are expensive to solve. Minimisation proce- 
dures can also lead to a more stable formulation than 
elliptic formulations, but arc usually slower. I will re- 
strict myself to the solution of the elliptic equation. As 
described below, I have also tried a simple Jacobi solver 
that is conceptually identical to a flow method, and it 
shows the expected robustness. 

It is possible to change the nature of equation (1) by 
extending it into the three-dimensional time slice. This is 
e.g. done by changing from the unknown function h(9, </>) 
to a level set h\(9, </>), where A is the level set parameter 
[13]. The advantage of such an extension is that multi- 
ple horizons can be found at the same time, and that no 
initial guess for the horizon location is needed. The dis- 
advantage is that such an extension will lead to a slower 
implementation because of the additional dimension. 



III. COORDINATES 

The above choice to use spherical coordinates to pa- 
rameterise the horizon introduces a preferred coordinate 
system into the otherwise coordinate-independent equa- 
tions (1), (2), and (3). Another preferred coordinate sys- 
tem comes from the grid used in the time evolution code 
for the whole spacetime. Often, gij and Kij are given in 
Cartesian components on a Cartesian grid. This raises 
the question as to which coordinate system to use in the 
numerical implementation of the horizon finder. The co- 
ordinate choice will of course not influence the physical 
results, but it can make the implementation more com- 
plicated if it requires interpolation, or less accurate if it 
leads to coordinate singularities. It can also make the 
results more or less difficult to interpret, as tensor quan- 
tities are numerically always expressed through their co- 
ordinate components. Tensors in different coordinate sys- 
tems cannot easily be compared, even if they are given 
at the same grid points. 

I chose to represent the quantities on the horizon us- 
ing their Cartesian tensor components, even when they 
are given on grid points on the horizon surface. This has 
the advantage that the tensors <?y and do not have 
to have their components transformed, although they do 
need to be interpolated to the apparent horizon location. 
It also has a disadvantage, because the partial deriva- 
tives that are calculated on the horizon are given only 
in spherical components, and those have then to be ex- 
plicitely transformed into Cartesian components. 

The overall advantage of using Cartesian tensor com- 
ponents is that non-scalar quantities on the horizon can 
more easily be compared to quantities given in the whole 
spacetime, and that the coordinate singularities at the 
poles do not influence the representation of such quan- 
tities. The singularities will of course still influence the 
calculation of quantities on the horizon, but the represen- 
tation of the results will not have coordinate singularities 
any more. 
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Figure 1: North pole of the grid on the horizon, in a polar 
projection. This shows the grid lines as they lie on the sphere. 



Let me use the indices i, j, k for Cartesian components 
(i, j, k e [x, y, z]), u, v, w for spherical components in 3D 
(u,v,w € [r, 6, <fi]), and a, b, c for spherical components 
on the horizon (a,b,c e [9,(j>]). This convention distin- 
guishes clearly which quantities are actually defined on 
what manifold, and also makes coordinate transforma- 
tions explicit. 



IV. DISCRETISATION 

The discretisation scheme of the 3D quantities, i.e. 
those living in the whole time slice, does not play a role 
when locating a horizon. There only has to be a way 
to interpolate these quantities onto the horizon surface. 
These quantities are usually discretised on a Cartesian 
grid. It is also possible to use spherical coordinates, or 
to use mesh refinement, without influencing the way in 
which the apparent horizon function H is evaluated. 

I chose to discretise the 2D quantities, i.e. those quan- 
tities living on the horizon, by using a polar #-0-grid with 
constant grid spacings d6 and d(j). A constant grid spac- 
ing works well when the horizon has a shape that is not 
too far from a sphere. Experiments show that distorted 
(peanut-shaped) horizons still work fine. 

Polar coordinates create coordinate singularities at 
9 = and 9 = 7r, which I avoid by having the grid 
points staggered with respect to the poles. I find that 
putting scalar quantities on such a grid generally works 
fine and requires no special treatment near the poles. 
On the other hand, naively putting tensor components 
in spherical coordinates on such a grid is a bad idea, as 
they either become zero or diverge at the poles. 

The boundary condition in ^-direction is periodicity: 
f(9, </>) = f(9, (f> + 2ir) for any function / living on the 
horizon. The boundary condition in ^-direction, i.e. 
across the poles, is a bit more involved. It is /(#, (j)) = 
P ■ f(—9, 4>-\-tt) for arbitrary tensor components /, where 
the parity P = (—1)** depends on the rank r of the tensor 
of which / is a component. Figures 1 and 2 demonstrate 
the polar boundary condition, especially how the shift by 
7r in ^-direction comes about. 
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Figure 2: North pole of the grid on the horizon, as seen from 
the grid. One can see that the black and white dots have 
multiple images on the coordinate grid, separated by ir in the 
^-direction. 



This grid allows numerical partial derivatives to be 
taken only in the 8- and (^-directions. The apparent hori- 
zon equations (1), (2), and (3) from above also need r- 
derivatives, and those cannot be taken numerically. How- 
ever, the r-derivative of the coordinate transformation 
operators are known analytically, and one can see from 
the definition of F that d r F = 1 by construction. Fi- 
nally, the partial derivative of the three-metric, which is 
given in Cartesian coordinates, is calculated before the 
metric is interpolated onto the horizon. 

My choice of evaluating the apparent horizon function 
directly on the grid that forms the horizon surface seems 
to be an uncommon one. It is also possible to evaluate 
the apparent horizon function not on the surface itself, 
but instead in the three-dimensional spacetime on those 
grid points that are close to the horizon [11, 14, 16]. This 
method requires interpolating in both directions between 
the surface and the time slice. It has the advantage that 
it is basically independent of the apparent horizon dis- 
cretisation, but it is also more expensive. In addition, the 
domain of the equation, i.e. the set of active grid points, 
changes with the horizon surface, leading to further com- 
plications. 

While I represent the apparent horizon shape using an 
explicit grid, one different commonly used way is to ex- 
pand the horizon surface location in spherical harmonics 
[1, 2, 8, 10]. Compared to using a grid on the hori- 
zon surface, an expansion in spherical harmonics has the 
advantage of having no coordinate singularities. Fur- 
thermore, one has direct insight into and control over 
the high spatial frequency components. Such an expan- 
sion would lend itself naturally to a multigrid algorithm. 
A multipole representation is by principle restricted to 
star-shaped regions, while an explicit representation is 
generic; however, this disadvantage is not important in 
practice. More problematic is that integrations over the 
surface are rather expensive, as they scale with 0(^ ax ), 
were / max is the number of multipole moments used. 
With an explicit grid, these integrations scale with 0(n), 
where n is the number of grid points, and where one 
should choose n — O(^ax) f° r a f an " comparison. 



V. EVALUATING THE APPARENT HORIZON 
FUNCTION H 

The apparent horizon finder uses the algorithm de- 
scribed below to evaluate the apparent horizon func- 
tion H. This description lists the important interme- 
diate quantities used in the apparent horizon finder, and 
explains the order in which they are calculated. This 
method to evaluate H is then subsequently used to find 
a solution for h. 

Although the quantities F, s l , and H are defined ev- 
erywhere in space, they only need to be evaluated on the 
horizon, i.e. on the surface determined by h. They still 
need to be defined in a neighbourhood of the horizon, 
so that their r-dcrivatives can be taken. Their continu- 
ation off the horizon is chosen arbitrarily, and is implic- 
itly determined by the above choice of spherical coordi- 
nates. Thus the finder calculates quantities on the two- 
dimensional horizon surface only, although the quantities 
live in three dimensions. 

It would be a bad idea to calculate the derivative of s l 
numerically, as s l is itself already calculated as a deriva- 
tive. Taking a derivative of a derivative leads to an effec- 
tive stencil that is larger, and hence slower, and that con- 
tains elements with a weight of zero. These zero weights 
lead to an odd-even decoupling of the grid points, which 
in turn leads to large errors. Instead, one uses the chain 
rule, and calculates VjS* directly from second derivatives 
of F. This makes it necessary to calculate numerically 
VjVj-F and ViS 1 along with V,F and s l . 

I describe explicitely which quantities are actually 
stored, and in which coordinate system the tensor quan- 
tities are calculated. 



Algorithm 

1. Give, as fixed background quantities, the three- 
metric gij{x k ) and the extrinsic curvature Kij{x k ) 
in the whole time slice. These quantities do not 
change while the horizon is located. 

2. Determine a 2-sphere h(x a ) for which the apparent 
horizon function is to be evaluated. This is done by 
an elliptic solver as described below, and is outside 
the scope of this algorithm. 

3. Calculate the horizon grid point locations in Carte- 
sian coordinates as x l (x a ), using the definition of 
h. 

4. Interpolate the three-metric, the extrinsic curva- 
ture, and the partial derivative of the three-metric 
onto the horizon. These quantities can conceptu- 
ally be viewed as representing gij{x u ) and Kij(x u ), 
but are stored as gij (x a ), Kij(x a ), and d k gij(x a ). 
(The partial derivative is needed later.) Note that 
these tensors now live on the horizon, but have their 
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components still in Cartesian components, prevent- 
ing coordinate singularities near the poles. 

5. Define the implicit horizon location F(x u ) = 
F(r, 9, (j>) = r — h(9, </>). This quantity is not stored 
explicitely. 

6. Calculate an outward-directed covector V v F(x u ) 
that is orthogonal to the surface. By definition, 
d r F — 1, and d a F = —d a h. This is stored as 
d v F(x a ) and d v d w F(x a ). (The second derivatives 
are needed later) . The calculation needs the values 
dbh(x a ) and dbd c h(x a ), which are calculated on the 
horizon surface. 

7. Dehnc the transformation operator T" that trans- 
forms from spherical to Cartesian covariant coordi- 
nate components. This operator is given by T" = 
dx u /dx l , where the x l {x u ) are given by the usual 

x = r sin(#) cos(0) 
y = r s'm(9) sin(</>) 
z = rcos(9) 

Evaluate this operator and its derivative at r = ft, 
and store it as T?(x a ) and djT^{x a ). 

8. Transform the components of the outward-directed 
covector into spherical coordinates: ViF(x u ) — 
TVW v F{x u ). Store W l F{x a ) and \7 t \7 3 F{x a ). 
The calculation of the second derivative needs the 
derivative of the transformation operator T, and 
also needs the partial derivative of the three-metric 
for the connection coefficients. 

9. Normalise the outward-directed covector ViF(x u ) 
and raise its index to s l (x u ), using the three- 
metric gij. Calculate its derivative ViS l (x u ). This 
needs the second derivative of F. Store s l (x a ) and 
V,s*(x a ). 

10. Calculate the apparent horizon function H = 
V t s l + K tJ (s l si - gv). Store it as H(x a ). 

For the above calculation I use a second-order uniform 
Cartesian interpolator, and second-order centred differ- 
ences to approximate the partial derivatives. Altogether, 
the apparent horizon function H thus depends on the 
background spacetime as given by and Kij, and de- 
pends on the horizon location ft. 

VI. FINDING THE HORIZON 

In order to solve a nonlinear equation, one needs initial 
data, which is in this case an initial guess for the horizon 
location. This initial data selects which apparent horizon 
will be found, if there are several present in the time slice. 
During the evolution one can easily use the location from 
the previous time slice as initial data; this works very well 



in practice and is called apparent horizon tracking. How- 
ever, finding horizons in an unknown spacetime is more 
difficult. Depending on how e.g. the initial configuration 
for an evolution run is constructed, one does not know 
whether, where, or of what shape the horizons are. A 
similar problem is encountered during the evolution, as 
new apparent horizons may appear. One does not know 
their location and not even their time of appearance in 
advance. 

Shoemaker ct al. [14] have implemented a mechanism 
in their apparent horizon finder that automatically de- 
tects whether a two black hole system has a common 
horizon, and if not, finds the two separate horizons in- 
stead. I did not use such a mechanism; if a horizon does 
not exist, I had my solver abort. I found it to be sufficient 
in practice to set up the finder manually with different 
initial guesses, corresponding to the different suspected 
horizons in a time slice. 



A. Newton-like method 

I implemented a fast Newton-like method by calling the 
PETSc library [6, 7]. This library provides several effi- 
cient parallel solvers for nonlinear elliptic equations. It 
uses Newton-like methods to reduce the nonlinear prob- 
lem to a linear one, and then offers Krylov subspace 
methods to solve it. 

Additionally, PETSc offers to calculate numerically the 
Jacobian that is necessary for a Newton-like method. I 
estimate that this is probably about an order of magni- 
tude slower than a hand-coded Jacobian; it is therefore 
not recommended by the PETSc authors. On the other 
hand is it much easier to implement, and most impor- 
tantly also much easier to get correct. Thornburg [16] 
explains the necessary and tedious steps of implementing 
an explicit Jacobian in detail. Huq et al. [11] use numer- 
ical perturbations instead, and their explicit calculation 
of those seems to complicate their implementation signif- 
icantly. Fortunately, PETSc is fast enough in my case, 
so I never tried a hand-coded Jacobian. Tracking a hori- 
zon with a reasonable resolution takes only about half a 
second on a notebook with yesteryear's hardware. 



B. Flow-like method 

Unfortunately, the apparent horizon equation is rather 
nonlinear, and the radius of convergence of Newton- 
like methods seems to be very small. Instead of a fast 
Newton-like method, one can apply a slow but robust 
flow-like method, which can be viewed as essentially an 
iterative Jacobi-type relaxation method. Such a method 
has a larger radius of convergence, and is in practice reli- 
able in finding the outermost apparent horizon. For this 
method one generally uses a large sphere as initial data. 
The horizon shape then flows gradually inwards, follow- 
ing the expansion of the horizon surface [14]. 
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With the above coordinate choice, the grid spacing 
near the poles is very small, and this requires very small 
convergence factors of the order of 10~ 3 or less for a rea- 
sonable resolution. Finding a horizon with a coarse res- 
olution takes several minutes on the above hardware. In 
principle one could first flow-find a coarser approxima- 
tion to the horizon, and then Newton-find the horizon 
with a higher resolution, but I have not tried this. 

I do not think that a multigrid algorithm can be used 
to accelerate flow-finding, because the coarse resolution 
is already so coarse that no coarser grid could reasonably 
describe a spherical shape. It might be worthwhile to use 
a different discretisation scheme for the horizon surface, 
so that the grid spacing does not tend to zero near the 
poles, or to use spherical harmonics as mentioned above. 



run 


dx d<p 


R/2 


Madm 


a 


1/4 2tt/36 


0.991111 


1.00043 


b 


1/8 2tt/72 


0.997693 


1.00010 


c 


1/16 2tt/144 


0.999429 


1.00002 







1.0 


1.0 



Table I: Horizon and ADM masses for a black hole with M = 1 
and a = 0. Both the spacetime and the horizon surface vary 
in resolution. The last row gives the analytic values. 



runs 


R/2 Madm 


a, b, c 


3.79 4.43 


a, b 


3.85 4.37 


b, c 


4.04 4.18 



Table II: Convergence factors for the results from table I. A 
value of 4 indicates perfect second order convergence. 



VII. ANALYSING THE HORIZON 

Once the location of the horizon is known, one can 
calculate its area A by integrating over the surface, and 
from that the horizon radius R via the definition A = 
AttR 2 . The irreducible mass is given by R/2, which is 
different from the total mass M in case the spin a is 
non-zero. 

One can also calculate the lengths of the equatorial 
and two polar circumferences (at (f> — and at <f> = tt/2). 
These circumferences can give a hint as to the coordinate 
distortion of the horizon. Of course, these hints are only 
reliable if one knows in advance about the symmetries of 
the horizon, which is only the case in testing setups. 

If one assumes that the horizon has a spin that is 
aligned with the z axis, then one can use the area and 
the equatorial circumference to calculate the spin mag- 
nitude. In the elliptic coordinates that one prefers for 
Kerr black holes (see e.g. [9, section 3.3]), a horizon with 
mass M an d spin a is located at the coordinate radius 
r = M + y/M 2 - a 2 . It has the area A = 47r(r 2 + a 2 ) 
and the equatorial circumference L = 2n(r 2 + a 2 )/r. (I 
thank Badri Krishnan for pointing out the latter to me.) 
These equations can be solved for a and M. 

It turns out that the estimate for the total mass M that 
is calculated this way has a reasonable accuracy, while 
the estimate for the spin a is rather inaccurate when a is 
small. In this case, a 2 can even become negative. This is 
mainly caused by accumulation of numerical errors. An 
O(10 -2 ) error in a 2 shows naturally up as an O(10 _1 ) 
error in a. Of course, the values for a and M are only 
correct when the spin is indeed aligned with the z axis, 
and that will usually not be the case (nor will it be easily 
verifiable) in a spacetime that is given numerically. How- 
ever, if the horizon should be isolated, then the isolated 
horizon formalism [3, 4] can provide a way of reliably cal- 
culating the spin and thus the total mass of the horizon. 



VIII. EXAMPLES 

I have implemented this apparent horizon finder al- 
gorithm within the Cactus framework [5], which has re- 
cently gained some popularity in the numerical relativity 
community. This will potentially make it easier for other 
people to use my implementation of this algorithm not 
only as an external programme, but also directly from 
within their code. 

Below I give results from applying this finder to con- 
figurations with one and multiple horizons. First I use 
single black holes in Kerr-Schild coordinates as test data 
to judge the accuracy and quality of this method. Then 
I apply the finder to superposed Kerr-Schild data to 
demonstrate the behaviour of this method in the pres- 
ence of multiple black holes, and compare the horizon 
masses to the corresponding ADM quantities. 

A. Single black hole 

I performed a convergence test with Kerr-Schild data. 
(See e.g. [9, section 3.3.1] for a definition of the Kerr- 
Schild coordinates.) Table I contains results from runs 
with a black hole of mass M = 1 and spin a = 0. The 
resolution of both the 3D time slice and the 2D grid on 
the horizon increases downwards. The last line contains 
the analytic values. R/2 is the irreducible horizon mass, 
and Madm is the ADM mass (see e.g. [17, eqn. (7.6.22)]) 
of the whole spacetime, which I calculated numerically 
for reference. 

In general, I would consider the resolutions 1/4, 1/8, 
and 1/16 to be coarse, reasonable, and fine, respectively. 
Corresponding time evolution simulations need to run in 
the year 2002 on a notebook, a workstation, or a super- 
computer. 

Table II shows the convergence factors for three-way 
and two-way convergence tests from the above runs. The 
convergence factors / are calculated in the usual way via 
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run 


d<t> 


R/2 


a z 


M 


A 


2tt/36 


0.960886 


0.496507 


0.994655 


B 


2tt/72 


0.964675 


0.498183 


0.998537 


C 


27T/144 


0.965594 


0.498815 


0.999511 


D 


27r/288 


0.965819 


0.498997 


0.999754 







0.965926 


0.5 


1.0 



Table III: Horizon masses and spins and ADM masses for a 
black hole with M = 1 and a z = 1/2. Only the horizon sur- 
face resolution varies. The last row gives the analytic values. 



A, B, C 

B, C, D 



R/2 a z 



M 



4.13 2.65 3.99 
4.08 3.47 4.01 



Table IV: Convergence factors for the results from table III. 
A value of 4 indicates perfect second order convergence. 



/ = (C - M)/{M - F) or / = (C - A)/{F - A), where 
C, M, F, and A represent the coarse, medium, fine grid, 
and analytic values, respectively. A factor of 4 indicates 
perfect second-order convergence. The results indicate 
that a spatial resolution of dx = 1/4 is not sufficient to 
be in the convergence regime, and also that the numerical 
ADM masses reported here might not be very accurate. 

In another convergence test, shown in table III, I 
kept the spatial resolution of the time slice constant at 
dx = 1/8, and varied only the resolution of the apparent 
horizon grid. This time the black hole had a mass M = 1 
and a spin a z — 1/2, again in Kerr-Schild coordinates. 
The apparent horizon spin a z is estimated via its equa- 
torial circumference. The last line of the table shows the 
analytic values. 

The convergence factors for this test are shown in table 
IV. As the spacetime is now fixed, and differs slightly 
from the analytic solution due to the discretisation errors, 
I omit the two-way convergence tests. It is clearly visible 
that the spin estimate does not converge to second order. 
This is unfortunate, but not of much consequence, as 
this method for the spin calculation is not applicable in 
general anyway. 

Figure 3 shows horizon shapes as reported by the 
finder. It shows apparent horizons of four black holes 
with varying spins and boosts. The singularity is at the 
origin for a = 0, and is a ring with radius a z in the 
x-y plane when the black hole is spinning. This figure 
shows nicely that excision boundary conditions become 
more complicated for spinning black holes, as the appar- 
ent horizon and the singularity can get rather close to 
each other. 



B. Multiple black holes 

In order to prepare and analyse initial data for binary 
black hole collision runs, I superposed two Kerr-Schild 
black holes as proposed in [12], and then solved the con- 
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Figure 3: Apparent horizons for Kerr-Schild black holes with 
varying spins and boosts. Displayed is a cut through the x-z 
plane. 



straint equations. This results in spacetimes with either 
two separate or a single distorted ("merged") black holes. 
One can freely specify the locations, masses, and spins of 
the black holes prior to the superposition. The relation 
to the properties of the superposed black holes is not 
known, but it is hoped the superposition will not change 
them much. 

For the superposition, I chose a grid spacing of dx — 
1/4, an outer boundary that has a distance 8.0 from the 
centres of the black holes, and excised a spherical region 
with radius 1.0 around the singularities. I superposed the 
three-metric and the extrinsic curvature without atten- 
uation. I then solved the constraint equations using the 
York-Lichnerowicz method with a conformal transverse- 
traceless decomposition. (This method is described e.g. 
in [9], section 2.2.1.) I used Dirichlet boundary condi- 
tions for the vector potential and the conformal factor, 
keeping the boundary values from the superposition. 

Whether or not a common apparent horizon exists is a 
hint as to whether the black hole configuration contains 
results from two separate or a single merged black holes; 
this fact is not known otherwise. The only way to find out 
would be to evolve this configuration for a long enough 
time so that the event horizon could be tracked back- 
wards in time. This is unfortunately not easily possible 
(if at all) with today's black hole evolution methods. 

I created a series of initial data configurations for two 
equal-mass (M = 1) non-spinning non-boosted black 
holes with varying distances. Table V shows the com- 
mon horizons masses Mq = Rc/2, the individual hori- 
zons masses Mj = Ri/2, and the ADM masses of the 
spacetimes Madm- The spins of the superposed black 
holes are known to be zero for symmetry reasons. 

The row z — 0.0 contains two black holes that were 
superposed at the same location, i.e. form a single black 
hole in a coordinate system different from Kerr-Schild. 
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z 


Rc/2 


Ri/2 Madm 


E B 


Er 


0.0 


2.280 


— 


2.306 





+0.026 


1.0 


2.245 


? 


1.973 


? 


-0.272 


1.75 


2.167 


1.137 


1.986 


-0.107 


-0.181 


2.0 


2.130 


1.087 


1.989 


-0.044 


-0.141 


2.5 


2.023 


1.041 


1.994 


-0.059 


-0.029 


3.0 


(2.101) 


1.020 


1.997 


(+0.060) (-0.104) 


4.0 




1.002 


2.000 




0.004 


6.0 




0.991 


2.000 


+0.018 


8.0 




0.989 


2.000 


+0.022 



Table V: Masses of the common horizons, the individual hori- 
zons, and the ADM masses for two superposed black holes 
for varying distances. Eb is binding energy of the two black 
holes, and Er is amount of radiation present in the spacetime, 
as explained in the text. For z = ±3.0, no common horizon 
was detected any more. The last three rows show the sum of 
Eb and Er only. 
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Figure 4: Common and individual horizons for superposed 
Kerr-Schild data at z — ±2.0. 



This is a consequence of the way the extrinsic curvatures 
are superposed. Note that the mass is not simply twice 
the mass of a single black hole. Also, as this spacetime 
is not spherically symmetric due to the outer boundary 
condition, the ADM mass and the horizon mass need not 
be the same. 

In the row z = ±3.0, the apparent horizon finder did 
not converge, but was close. The L 2 - norm of the residual 
was still 0.0248 when the solver aborted. That means 
that from this row on there is no common horizon any 
more, but the solution error here is still comparable to 
the discretisation error. 

The difference Eb — Mc — 2Mi can be interpreted as 
the binding energy of the system. (Here Mc — Rc/2 
and Mi = Rj/2 because the spin is zero.) The binding 
energy magnitude seems to decrease with distance, and 
is negative as it should be, except in the row z — ±3.0 
where there was no common horizon detected any more. 

The difference En — Madm — Mc is the amount of 
radiation present in the spacetime. Of course, En must 
be positive, meaning that the negative values found nu- 
merically are unphysical and point to numerical errors. I 
suspect that the error lies mainly with the calculation of 
the ADM mass. I calculate the ADM mass as a surface 
integral near the outer boundary. When I evaluate the in- 
tegral at locations further inwards, closer to the apparent 
horizon, then the ADM mass increases and gets closer to 
the horizon mass. This is inconsistent and cannot happen 
when the constraints are satisfied. However, in a numer- 
ically given spacetime the constraints are only satisfied 
up to the discretisation error. It is thus in my opinion a 
bad idea to calculate the ADM mass (which is a global 
quantity) from information near the outer boundary only, 
relying on the constraints to propagate the necessary in- 
formation to the boundary. A better method to calculate 
the ADM mass is needed. I decided to present the ADM 
values here nonetheless, because they (although inaccu- 
rate) still provide some insight into the spacetimes. 
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Figure 5: Common and individual horizons for superposed 
Kerr-Schild data at z — ±1.75. The individual horizons do 
overlap, and are surprisingly only minimally distorted. 

Figure 4 shows the apparent horizons of two super- 
posed black holes with centres at z = ±2.0. The common 
apparent horizon indicates that there is a common event 
horizon present, and that this time slice thus contains a 
single, merged black hole. In figure 5, the black holes are 
closer together at z — ±1.75. The two inner apparent 
horizons actually do overlap; this is not a numerical er- 
ror. It is surprising to see that the shape of the horizons 
is not distorted more. But unlike event horizons, there is 
no reason why apparent horizons should deform or join 
when two black holes approach each other. 

Figure 6 shows the the common apparent horizons for 
a series of black holes with increasing distances. As the 
black holes are placed further apart, the common hori- 
zons become unsurprisingly more elongated. Figure 7 
shows both the individual and common black holes for 
varying distances for comparison. Near z — ±3.0, the 
common horizon ceases to exist. 
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Figure 6: Common horizons for a series of merged black holes 
with increasing distances. 



IX. SUMMARY 



Apparent horizons provide valuable insight into space- 
times, such as mass and spin estimates for black holes. 
They are indicators for event horizons, and help keep 
track of the location of singularities. They provide nec- 
essary information for the time evolution of black hole 
spacetimes. This apparent horizon finding algorithm pro- 
vides a fast way to track horizons during a numerical evo- 
lution, and also a robust way to find horizons in initial 
data, or as they appear during an evolution. 
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Figure 7: Individual and common horizons for superposed 
Kerr-Schild black holes near z = ±3.0, where the common 
apparent horizon disappears. 
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